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SUMMARY 

A  Fortran  program  has  been  developed  for  the  Maximum  Likelihood  Estimation  of 
parameters  in  non-linear  systems.  The  program  structure  uses  subroutines  to  describe 
the  problem  and  define  problem-specific  elements,  while  the  main  program  is  designed  to 
be  problem-independent  as  far  as  possible.  In  the  present  note  an  application  of  the  program 
to  compatibility  checking  of  aircraft  dynamic  flight  test  data  has  been  studied.  Using 
simulated  time  histories  of  longitudinal  manoeuvres,  the  conditions  for  satisfactory 
performance  have  been  identified.  It  has  been  shown  that  for  a  practical  manoeuvre  shape, 
record  length  and  sampling  rate  and  for  reasonable  noise  levels  on  the  measured  data, 
instrument  errors  can  be  identified  to  good  accuracy  and  a  set  of  compatible  time  histories 
reconstructed. 
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1.  INTRODUCTION 

The  Aircraft  Behaviour  Studies — Fixed  Wing  Group  has  been  actively  involved  in  the 
determination  of  aircraft  aerodynamic  characteristics  from  flight  test  data  for  a  number  of  years. 
Good  results  have  been  obtained  (e.g.  Ref.  1)  when  the  mathematical  model  describing  the 
aircraft  motion  is  reasonably  well  known  and  linear  or  almost  linear.  In  this  case  powerful 
methods  are  available  (e.g.  Ref.  2)  which  can  identify  the  aerodynamic  parameters  in  the  presence 
of  both  state  (or  process)  and  measurement  noise  including  instrumentation  biases. 

If  the  model  is  essentially  non-linear  then  the  problem  becomes  more  difficult  and  is  further 
complicated  by  the  fact  that  such  models  are  often  not  well  defined.  In  particular  this  applies  to 
the  dynamics  of  aircraft  at  high  incidence.  Reference  3  provides  a  review  of  high  angle  of  attack 
aerodynamics  and  points  to  the  substantial  effort  in  this  area  over  recent  years.  Nevertheless 
the  present  state  of  knowledge  is  still  largely  descriptive.  At  the  same  time  methods  for  identi¬ 
fication  of  aerodynamics  in  high  angle-of-attack  regimes  are  in  a  relatively  early  stage  of  develop¬ 
ment  and  no  standard  approach  has  been  established.  Reference  4  is  an  example  of  one  approach 
which  makes  clear  the  difficulty  of  the  problem  due  to  the  presence  of  noise  and  instrument 
errors  in  the  measurements,  and  uncertainties  in  the  (non-linear)  model  structure.  Because  of 
this  it  would  be  desirable  to  separate  out  these  two  sources  of  error.  This  can  conveniently  be 
done  in  the  aircraft  case  by  separating  the  kinematics  from  the  full  dynamical  description  of  the 
motion.  The  kinematic  equations,  relating  accelerations,  velocities  and  displacements  are  non¬ 
linear  but  well  defined  and,  ideally,  can  be  analysed  to  provide  the  unknown  instrumentation 
biases  and  an  optimal  estimate  of  the  state  and  measurement  vectors.  The  resulting  “error-free” 
records  are  then  used  in  a  full,  but  not  necessarily  well-defined,  dynamic  model.  The  methods  of 
regression  analysis,  for  example,  can  conveniently  be  applied  in  attempting  to  identify  the 
parameters  as  well  as  the  structure  of  this  model.  Regression  analysis  can  be  expected  to  give 
good  results  when  measurement  bias  errors  are  absent  (Ref.  8).  This  two-stage  approach  can, 
of  course,  be  applied  equally  well  to  the  linear  well-defined  model  case  as  outlined  in  References 
6,  7  and  8.  In  those  references  an  Extended  Kalman  Filter  is  used  in  the  first  stage,  i.e.  for 
estimation  of  the  aircraft  dynamic  state  and  of  instrument  systematic  errors  via  the  aircraft 
kinematic  equations.  Reference  6  obtains  good  results  with  a  restricted  number  of  variables 
while  relying  on  high  quality  instrumentation.  References  7  and  8  present  a  broader  formulation 
but  do  not  discuss  the  requirement  for  successful  implementation  such  as  noise  levels  and  the 
need  to  specify  them,  record  length,  sampling  rates,  etc.  Reference  9  points  to  some  doubt  on 
the  suitability  of  an  Extended  Kalman  Filter  as  a  parameter  estimator  (for  identifying  systematic 
instrument  errors),  then  proceeds  to  apply  a  maximum  likelihood  method,  ignoring  process 
noise,  and  presents  some  preliminary  results.  The  neglect  of  process  noise  in  the  kinematic 
equations  can  be  justified  on  the  basis  that  it  is  small  since  instrumentation  used  for  the  measure¬ 
ment  of  the  inputs  is  the  most  accurate  available.  Furthermore,  modelling  errors  are  absent 
since  the  kinematic  equations  can  be  considered  to  be  exact.  Reference  10  notes,  for  their  restricted 
case  and  accurate  instrumentation,  that  application  of  the  Maximum  Likelihood  algorithm  (no 
process  noise)  and  the  extended  Kalman  Filter/Smoother  have  been  found  to  yield  similar  results. 

In  this  note  a  Maximum  Likelihood  (ML)  program  is  developed  for  application  to  a  general 
non-linear  system  with  no  process  noise.  The  user  is  required  to  specify  the  system  under  study 
via  two  subroutines,  one  setting  out  the  differential  equations  of  the  system  and  the  other  speci¬ 
fying  the  system  sensitivity  matrix.  Other  problem-dependent  subroutines  are  used  for  initialisa¬ 
tion  and  for  calculation  of  the  output  responses.  The  program  has  been  applied  here  to  the 
estimation  of  states  and  instrument  systematic  errors  using  the  aircraft  kinematic  equations, 
as  outlined  above.  The  aim  of  this  study  is  to  assess  the  performance  of  the  ML  algorithm  for 
this  problem,  to  establish  the  conditions  required  for  its  successful  implementation  and  to 
provide  a  baseline  for  comparison  with  other  methods.  To  this  end  the  effects  of  noise  levels, 
record  lengths  and  sampling  rates,  manoeuvre  shape  and  the  number  of  time  histories  matched, 
etc.,  have  been  studied  and  the  results  presented. 
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A  theoretical  description  of  the  method  including  an  outline  of  the  ML  program  is  presented 
in  Section  2  and  results  of  its  application  to  the  aircraft  kinematic  equations,  using  simulated 
data,  are  given  in  Section  3. 


2.  DESCRIPTION  OF  METHOD 

A  brief  theoretical  discussion  of  the  Maximum  Likelihood  method  is  followed  by  a  more 
detailed  discussion  of  its  application  to  the  aircraft  kinematic  equations.  Consideration  of  pos¬ 
sible  strong  correlations  among  some  of  the  parameters  leads  to  several  possible  alternative 
approaches.  A  description  of  the  computer  program  ends  this  section. 

2.1  Theory 

The  system  is  assumed  to  be  described  by  a  set  of  non-linear  dynamic  equations  of  the  form 

x(t)  =  f(x(r),u(/),£),  (1) 

z(/i)  =  g(x(/i).u(f<).£)+n('<).  (2) 

where  x  is  the  state  vector, 
u  is  the  input  vector, 
z  is  the  observation  vector, 
n  is  the  measurement  noise  vector, 

£  is  the  vector  of  unknown  parameters. 

The  state  equation  (1)  is  assumed  to  be  free  from  process  noise  and  the  measurement  noise 
is  assumed  to  be  zero  mean  white  Gaussian  noise.  The  measurements,  z,  are  made  at  a  finite 
number  of  time  points.  It. 

In  order  to  estimate  the  unknown  parameters  a  likelihood  function,  L( z|f),  is  defined  which 
expresses  the  probability  of  obtaining  the  measurements,  z,  given  the  parameters,  £.  The  maxi¬ 
mum  likelihood  (ML)  estimate  of  £  is  the  value  which  maximises  L(z\£),  evaluated  with  the 
measured  responses.  The  ML  estimate  has  a  number  of  desirable  properties,  including  consistency 
and  efficiency  (see  Ref.  11).  Using  the  Gaussian  assumption  for  the  measurement  noise  the 
likelihood  function  for  the  present  case  can  be  written  (e.g.  Ref.  11): 

Hz |f)  =  [2ir*|R|  ]  "/2exp[-£  I  (z(/<) -z£(/())tR' '(*('<)-*£('<))]  (3) 

i-t 

where  m  is  the  length  of  the  observation  vector,  z, 

N  is  the  number  of  time  points. 

In  equation  (3)  z{  is  calculated,  for  a  particular  £,  from  equations  (1)  and  (2)  neglecting  the 
measurement  noise,  n,  while  R  is  the  covariance  of  the  residuals,  z(f<)— z c(/<). 

Taking  the  log  of  equation  (3)  for  simplicity,  the  maximisation  of  L  is  equivalent  to  the 
minimisation  of  the  cost  functional 

J(£, R)  =  JwNiog2^+JAlog]R|  +  J  Z(z(li)-zc(ti))TR-'{z(ti)-z((U))  (4) 

1 

For  given  R  the  first  two  terms  of  J  are  constant  and  the  following  simple  cost  functional 
results : 

-'(f)  =  i  I  (z(/<)  — z{(/())TR  '(zOO-zJtO)  (5) 

I  I 

an  estimate  for  R  can  be  obtained  by  minimising  J(£, R)  (eqn  (4))  with  respect  to  R.  This  results  in 

R  =  ‘  I  {z{U)-zc{t())(z(li)-Z({U))T.  (6) 

j>  | 

Thus  the  problem  splits  into  two  parts.  For  fixed  £  equation  (6)  maximises  the  likelihood 
function  with  respect  to  R,  while  for  given  R  the  cost  functional  given  by  equation  (5)  is  mini- 
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mised  to  estimate  j.  The  optimisation  algorithm  used  to  achieve  this  is  the  modified  Newton- 
Raphson  algorithm  as  documented  in  Reference  12.  Starting  from  an  initial  estimate  for  £, 
revised  estimates  are  obtained  iteratively  from 

=  6-[VA6)]-I[V(6)]r.  (7) 

The  first  gradient  matrix  is  obtained  directly  from  equation  (5): 

m  =  Z  (z(/<) -*f(/i))TR  1 F £(ze(/«))]  (8) 

i=i 

and  an  approximation  to  the  second  gradient  is  given  by 

V£W)  =  I  F£(z£(rO)]TR  ^£(z£(/<)).  (9) 

i=i 

The  term  V((zs(tt))  required  to  evaluate  equations  (8)  and  (9)  is  the  sensitivity  matrix  to  be 
obtained  from  the  system  equations. 

Thus  the  iterative  procedure  used  can  be  summarised  as  follows: 

(1)  start  with  initial  values  for  f,  R; 

(2)  with  fixed  R  use  equations  (7)— <9)  to  estimate  a  new  f ; 

(3)  from  equation  (6)  obtain  a  new  estimate  of  R  using  the  revised  £ ; 

(4)  repeat  steps  2  and  3  until  convergence. 

In  practice  a  few  iterations  are  done  with  R  fixed  before  step  3  is  included,  since  the  residual 
power  can  often  be  quite  large  in  the  first  iteration  or  two  and  hence  lead  to  a  worse  estimate  of  R 
than  the  starting  value. 

For  the  Gaussian  case  Reference  11  shows  that  the  maximum  likelihood  estimator  is 
asymptotically  consistent  and  efficient,  i.e.  for  large  N  the  estimates  for  f  are  normally  distributed 
about  the  true  value  with  covariance  given  by  the  Cramer-Rao  lower  bound.  Thus  a  measure  of 
the  accuracy  of  the  estimates  is  given  by  (Ref  11): 

covariance(f)  =  ^  ZJTj(z£(/f))]rR_1V£(z£(/j))  1  (10) 

The  right-hand  side  of  equation  (10)  can  be  seen  to  be  the  inverse  of  the  second  gradient 
matrix  as  approximated  by  equation  (9). 


2.2  Aircraft  Kinematic  Equations 

The  kinematic  equations  are  a  set  of  non-linear  equations  relating  the  position,  velocity  and 
acceleration  of  an  aircraft  with  reference  to  a  set  of  fiat  earth  axes.  The  system  can  be  defined 
as  follows: 

(i)  State  vector,  x  =  [u,ry\\<f>,6,ip,h]T 

The  state  equations  written  in  body  axes  (x,y,z)  fixed  in  the  aircraft  with  origin  at  the 
centre  of  gravity  can  be  considered  to  be  exact  (e.g.  Ref.  13): 

u  =  —  qw+rv+  ax— g  sin  0  ■>, 

v  =  —  ru+pw+ay+g  cos  9  sin  <f> 

w  =  qu—  pv+az+ g  cos  9  cos  <f> 

$  =  p-f<7sin  <f>  tan  #+rcos  <f>  tan  0  >  (11) 

6  =  q  cos  4>  -r  sin  <l> 

ijt  =  <7  sin  <£/cos  0-rr  cos  <£/cos  0 

h  =  u  sin  0—  »•  cos  0  sin  <f>  —  tv  cos  0cos  6 


where  u,  v,  w  are  linear  velocities  in  x,  y,  z  directions, 

<P,  6,  <p  are  roll,  pitch  and  yaw  attitudes 
h  is  altitude, 

p,  q,  r  are  roll,  pitch  and  yaw  angular  rates, 

ax,  ay,  az  are  linear  acceleration  in  x,  y,  z  directions, 

g  is  gravitational  constant. 

Note  that  the  h  and  p  equations  are  not  coupled  into  the  other  equations. 

(ii)  Input  vector,  u  =  [ax,ay,az,p,q,r]T 

In  general  measured  values  for  the  input  vector  are  corrupted  by  scale  errors,  instru¬ 
ment  bias  errors  and  random  noise.  For  example,  taking  the  pitch  rate,  q,  we  write 

q  =  (I  +  bq)qm-\-bq-\-nq,  (12) 

where  qm  is  measured  pitch  rate, 

A9  represents  scale  factor  error, 

bq  represents  instrument  bias  error, 

nq  represents  random  noise, 

with  similar  relations  for  the  other  elements  of  the  input  vector.  If  equation  (12)  and 
its  counterparts  are  substituted  into  equation  (11)  then  the  unknown  parameters  A„, 
bq,  etc.,  appear  explicitly  in  the  state  equation  and  in  addition  the  random  noise,  nq, 
gives  rise  to  state,  or  process,  noise.  However  accelerometers  and  gyros  used  in  measur¬ 
ing  the  input  quantities  are  the  most  accurate  instruments  used  and  random  noise 
levels  are  generally  small.  Thus  it  is  reasonable  to  assume  that  the  process  noise  is 
negligible. 

(iii)  Output  vector,  z  =  [V,pv,<xv,<p,0,ip,h]T 

Measurements  of  the  output  vector  are  corrupted  by  scale  factor  errors,  biases  and 
random  noise.  The  equations  for  calculating  the  outputs  are  taken  as: 

Fout  =  (I+Av)(u2-i-v2  +  w2)l'2  +  *v 

LM  «  J 

<£out  = 

^out  —  (1  -\-^e)Q-\-bg 
'pom  =  (1 +AJ>)</r-)-()^ 
bom  =  (1 

where  V  is  airspeed, 

«v.  |5v  are  incidence  and  sideslip  angles  at  the  sensor  (vane). 

The  term  <u  -pz  -)/u  in  the  (?v  equation  and  t'  e  equivalent  term  in  the  av  equation  are 
corrections  due  to  known  sensor  offsets  (■*>,», z„)  and  (xx,ya,Za)  relative  to  the  centre  of  gravity. 

Given  the  measurements  of  inputs  and  outputs  the  requirement  is  to  estimate  the  biases 
and  scale  factors  introduced  in  equations  (12)  and  (13).  In  addition,  although  the  state  equation 
is  noise  free,  the  initial  conditions  are  not  known  exactly  and  hence  become  part  of  the  unknown 
parameter  vector: 

f  =  [ baz ,  bay,  bat,  bp,  bq,  br,  bv,  bg,  bx,  b$,  bp, 

Aaz,  Aay,  Aa*,  Ap,  A9,  Ar,  Av,  A„  A^,  A^,  A^,  A*, 
u(0),  v(0),  h(0),  *(0).  0(O)F_ 


(14) 


The  initial  values  and  biases  for  yaw  attitude,  <f>( 0),  and  height  /r(0),  b/,  can  be  assumed 
to  be  zero  since  there  is  no  absolute  reference  implied  in  equation  (II).  The  number  of  para¬ 
meters  in  equation  (14)  can  be  reduced  if  some  of  the  biases  or  scale  factor  errors  can  be  assumed 
to  be  zero.  For  example,  References  7  and  8  neglect  all  scale  factor  errors  except  for  Av,  A^,  Xa 
while  Reference  9  retains  only  the  bias  errors. 

For  the  present  study,  the  problem  has  been  scaled  down  further  by  considering  the  subset 
for  longitudinal  motions  only  and  neglecting  all  scale  factor  errors.  The  reduced  system  can  be 
summarised  as  follows : 


(i)  Slate  vector,  x  =  [ u ,  w,  8,  h]T 

u  =  —(qm+b^w+iaxm+ba^—gsin  8 
W  =  (<7rn + bq)u +(azmJrbaz)-Vgco%  8 

&  =  (fm-\-bq 

h  =  u  sin  8—w  cos  8 


where  [axm,  azm,  <7m]T  is  the  measured  input  vector. 

(ii)  Output  vector,  z£  =  [F0„t,  *vout,  8 out,  AoutF 
Vout  =  («2  +  w2)>'2  +  /)v 


u  u 


8  out  =  8-\-bg 
h  out  =  b-\-bh 

* 

The  unknown  parameter  vector  to  be  estimated  is  now 

f  -  [bax,  baz,  V  bv,  b„  be,  u(0),  h(0),  0(O)F.  (16) 

The  initial  height,  /i(0)  and  the  bias  bn  can  be  assumed  to  be  zero  and  hence  are  not  included 
in  the  parameter  vector. 


2.3  Sensitivity  Matrix 

The  parameter  estimation  procedure  as  outlined  in  Section  2.1  requires  the  calculation  of 
the  sensitivity  matrix,  Ve(i ((tt)).  For  a  measurement  vector  of  dimension  m  and  parameter  vector 
of  dimension  p  the  sensitivity  matrix  is  of  order  m  by  p  with  the  (i,j)  element  being  the  partial 
derivative  of  the  rth  component  of  z£  with  respect  to  the  y'th  component  of  For  example  the 
(1,  2)  element  would  be  the  partial  derivative  of  Kout  with  respect  to  baz,  kout ,e>oi.  From  equation 
(15)  it  follows  that 

Fout.i >„  =  (m  .  Uba!+  W  .  H'6ar)/("2  +  M'2)1/2  ( 1 7) 

where  u,  w  and  the  partial  derivatives  and  are  all  functions  of  time.  The  velocity  com¬ 
ponents  u  and  >v  are  obtained  directly  from  the  integration  of  equation  (14).  At  the  same  time 
equations  for  maz  and  Wba!  can  readily  be  derived  from  the  state  equation  (14).  Thus,  taking 
partial  derivatives  of  equation  (14)  with  respect  to  baz  provides  the  required  set  (note  that 
V  =  0): 

Uba!  =  — (<jm~\'bq)Wl>aI  ) 

(18) 

7m-!-A)wt,„+l  J 

with  zero  initial  condition^.  .»  us  .  nations  (18)  have  to  be  integrated  simultaneously  with 
equations  (14)  to  enable  Vout,baz  to  be  evaluated  as  a  function  of  time  from  equation  (17). 

The  complete  sensitivity  matrix  for  the  present  system  is  summarised  in  Appendix  I  and  the 
full  system  of  differential  equations  requiring  to  be  integrated  in  order  to  calculate  the  elements 
of  the  sensitivity  matrix  is  presented  in  Appendix  2. 


5 


It  is  a  straightforward  process  to  extend  the  system  to  include  additional  parameters  if 
desired.  Alternatively,  if  change  of  altitude  is  not  measured  it  would  be  desirable  to  remove  h 
from  the  sytem  being  considered.  Thus  the  state  vector  and  output  vectors  would  reduce  to 
dimension  three.  The  sensitivity  matrix  would  reduce  from  4x9  to  3x9  and  the  total  number 
of  equations  to  be  integrated  would  reduce  from  the  23  in  Appendix  2  to  16. 


2.4  Correlated  Parameters 


Examination  of  the  output  equations  (15)  reveals  relationships  between  the  unknown 
initial  conditions  t/(0),  n(0)  and  0{ 0)  and  the  unknown  biases  bv,  b2  and  bH: 


V out(0)  =  («2(0)  f  »>  '(0))>  2  +  Z»v 
■(0)  (^m(0)  +  69) 

M(0)  ~  u( 0) 

0ou.(O)  =  0(0) +6,. 


«vou,(0)  =  tan 


ip'f 

L“0 


(19) 

(20) 
(21) 


For  initially  steady  level  flight  (<7(0)  =  qm(0)+bQ  =  0)  equation  (20)  is  slightly  simplified 
and  in  addition  a  relation  between  0(0),  bax  and  bai  follows  from  equations  (14): 


tan  0(0)  =  — 


tf.Vm(O)  +  bax 
fl-m(O)  +  baz 


(22) 


In  the  absence  of  measurement  noise  it  is  required  that  calculated  outputs  are  equal  to 
measured  outputs: 

Fout(0)  =  V m(0)  ) 

«vout(0)  =  «Vm(0)  (23) 

0out  =  0tn(O)  J 


Thus  any  change  in  initial  conditions  must  be  accompanied  by  a  change  in  one  or  more  of 
the  biases,  e.g., 

=  — A0(O),  etc.  (24) 

This  suggests  that  if  the  initial  conditions  are  estimated  then  the  related  biases  can  be 
calculated  directly  without  estimating  them  as  independent  parameters.  The  number  of  parameters 
to  be  estimated  is  thus  reduced.  However,  at  the  same  time,  to  be  consistent,  appropriate  elements 
of  the  sensitivity  matrix  have  to  be  modified.  For  example  equation  (24)  implies  that  bn  is  now 
a  function  of  0(0)  rather  than  an  independent  parameter  and  the  partial  derivative  of  b0  with 
respect  to  0(0)  is  —I.  Using  this  with  equation  (21)  gives  the  result  for  the  sensitivity  element 
0out,«(o)  =  0  rather  than  the  previous  value  of  1.  Similar  modifications  to  the  sensitivity  matrix 
result  when  bv  and  ba  are  taken  as  functions  of  the  initial  conditions  as  implied  by  equations 
(19)  and  (20).  Alternatively  it  is  possible  to  consider  u( 0),  u(0)  and  0(u)  to  be  functions  of  bv, 
bx  and  bH  rather  than  the  opposite,  thus  removing  the  initial  conditions  from  the  parameter 
vector  in  favour  of  bv,  h2  and  b„.  Finally,  for  initially  steady  level  flight  equation  (22)  suggests 
the  possibility  of  calculating  baz  (or  baz)  as  a  function  of  0(0)  and  baz  (or  baz ).  Modifications  to 
the  sensitivity  matrix  implied  by  these  options  are  summarised  in  Appendix  3. 

If,  in  the  noise  free  case,  all  the  parameters  are  estimated  as  though  fully  independent,  the 
correlations  between  initial  conditions  and  biases  as  expressed  in  equations  (19)  to  (22)  may 
well  be  expected  to  lead  to  less  accurate  or  biased  results.  On  the  other  hand,  in  the  presence 
of  measurement  noise  the  relations  given  in  equation  (23)  do  not  hold  precisely  so  that  the  pro¬ 
cedure  outlined  above  may  not  be  justified.  The  effect  of  noise  will  be  examined  further  in 
S  ’etion  3. 


2.5  Computer  Program 

The  basic  structure  of  the  program  is  shown  in  Figure  I.  In  the  first  stage  the  input  data 
provide  all  the  information  necessary  to  broadly  define  the  problem.  Thus,  for  example,  the 
number  of  equations  to  be  integrated  and  the  parameters  to  be  estimated  are  specified  as  well 
as  the  a  priori,  or  initial,  values  for  the  parameters  and  for  the  weighting  matrix  R.  In  addition, 
fixed  data,  c.g.  incidence  vane  location  (v,.  tv,,  r,),  etc.,  and  data  relating  to  computational  aspects 
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such  as  step  size,  total  number  of  time  points  and  maximum  number  of  iterations  to  be  performed 
are  provided  at  this  stage.  Finally,  time  histories  of  all  the  measured  inputs,  u(fi)>  and  outputs, 
z(/<),  are  read  in  and  stored  in  arrays,  since  they  will  be  repeatedly  called  upon  during  the  iteration 
process. 

The  iteration  loop  is  commenced  by  setting  up  the  arrays  and  matrices  to  be  used  in  the 
subsequent  calculation.  For  example,  arrays  representing  the  state  variables  and  their  derivatives 
with  respect  to  time  are  initialised,  together  with  the  sensitivity  matrix  and  the  first  gradient 
array  and  second  gradient  matrix  (equations  (8),  (9)).  The  array  and  matrix  formats  used  follow 
that  given  in  Reference  12  thus  enabling  many  of  the  utility  subroutines  provided  there  to  be 
used  in  the  current  program.  Much  of  the  matrix  manipulation  in  this  block  and  in  the  block 
labelled  “Estimate  New  Parameters”  closely  follows  that  documented  in  Reference  12.  For  the 
present  study,  in  order  to  allow  for  flexibility  in  changing  the  number  of  parameters  to  be  esti¬ 
mated,  and  hence  the  number  of  variables  to  be  integrated  (Appendix  2),  the  subroutine  1NIT 
has  been  written.  Given  the  parameters  to  be  estimated,  IN1T  works  out  which  equations  need 
to  be  integrated  and  sets  the  initial  conditions.  Finally,  before  entering  the  time  loop,  the 
integrator  is  initialised,  step  size  and  required  accuracy  set,  etc. 

The  time  loop  proceeds  to  integrate  the  necessary  equations  (  Appendix  2)  and  concurrently 
evaluates  the  sensitivities  (Appendix  1)  and  hence  the  gradients  (equations  (8),  (9))  as  well  as 
the  outputs  (equation  (15)),  all  based  on  the  current  parameter  values.  Integration  is  performed 
via  the  fourth  order  Runge-Kutta  routine  (subroutine  AN1NT)  although  a  predictor-corrector 
method  has  also  been  used.  The  integration  routine  calls  for  derivai.ve  values  which  are  supplied 
by  subroutine  DER1VS.  When  integration  to  the  next  time  point  is  complete,  sensitivities  and 
outputs  are  evaluated  via  subroutines  SENS  and  RESP  respectively  before  the  gradients  can  be 
calculated.  The  time  loop  proceeds  step  by  step  until  the  specified  number  of  points  is  achieved. 
Note  that  on  the  first  iteration,  measured  values  of  the  states  (or  values  deduced  from  the 
measurements)  are  used  in  calculating  the  gradient.  This  technique  is  useful  if  initial  parameter 
values  are  far  from  the  correct  values  (see  Ref.  12). 

All  the  information  required  for  the  estimation  and  updating  of  the  parameters  is  now 
available  (equation  (7))  and  this  is  now  done.  At  the  same  time  an  updated  weighting  matrix 
is  calculated  using  equation  (6).  As  explained  in  Section  2.1  this  step  is  only  carried  out  after 
the  first  few  iterations.  The  iteration  loop  proceeds  until  the  specified  number  of  iterations  is 
complete.  Alternatively,  it  would  be  possible  to  set  a  convergence  criterion  on  the  residuals 
but  this  has  not  been  done  here. 

Finally,  the  correlation  matrix,  whose  diagonal  elements  are  the  Cramer-Rao  bounds,  are 
calculated  using  equation  (10).  The  output  results  include  an  iteration  history  comprising  fit 
error  (residuals),  gradient,  determinant  of  second  gradient  matrix  and  new  values  of  estimated 
parameters  after  each  iteration.  The  final  values  of  the  parameters  and  their  Cramer-Rao 
bounds  are  then  summarised  followed  by  the  normalised  correlation  matrix.  A  second  output 
file  contains  time  histories  of  the  inputs  and  the  measured  and  calculated  outputs  to  be  used  for 
producing  time  history  plots  if  required. 

Since  non-linear  systems  can  vary  greatly  in  their  structure  it  is  difficult  to  write  a  general 
program  to  cater  for  all  possible  variations.  In  the  present  case,  although  the  main  interest  has 
been  in  the  aircraft  longitudinal  kinematic  equations,  an  attempt  has  been  made  to  separate  the 
problem-specific  elements  from  the  more  generally  applicable  program  elements.  Thus  the  sub¬ 
routines  INIT,  DERIVS,  SENS  and  RESP  are  problem-specific  and  would  have  to  be  altered 
for  each  different  system.  Some  alteration  in  the  problem  definition  box  (Fig.  1)  would  probably 
also  be  necessary.  Such  modifications  are  relatively  simple  and  make  it  possible  to  apply  the 
program  to  estimation  of  parameters  in  other  non-linear  systems  of  interest.  For  example,  it 
could  be  used  for  estimating  the  aerodynamic  parameters  in  a  non-linear  dynamic  model  of  the 
aircraft  using  the  “error-free”  responses  produced  by  the  present  program.  However,  for  that 
particular  problem,  regression  analysis  would  probably  be  more  convenient  to  apply. 

3.  RESULTS  AND  DISCUSSION 

The  longitudinal  system  defined  in  equations  ( 1 4)-(  1 6)  have  been  studied  in  some  detail 
using  simulated  data.  Measurement  time  histories  were  produced  by  an  auxiliary  program 
which  allowed  specified  biases  and  noise  levels  to  be  set.  Process  as  well  as  measurement  noise 
could  be  specified  in  order  to  study  the  effect  that  process  noise  had  on  the  results. 
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The  study  proceeds  systematically  by  first  looking  at  the  no  noise  case  and  assessing  the 
effects  of  changing  the  number  of  measured  time  histories  to  be  matched  and  the  performance 
of  the  different  procedures  (Section  2.4)  designed  to  cope  with  correlated  parameters.  This  is 
followed  by  results  with  measurement  noise  and  then  process  noise  added.  The  noise  is  in  all 
cases  white  and  Gaussian.  The  effects  of  record  length,  sampling  rate  and  noise  levels  are  investi¬ 
gated  and  finally  the  influence  which  manoeuvre  shape  has  on  the  results  is  briefly  examined. 
Some  care  has  been  taken  to  ensure  that  integration  errors,  due  to  too  large  an  integration  step 
size,  are  avoided  while  maintaining  step  size  as  large  as  possible  in  order  to  minimise  compu¬ 
tation  time.  Although  a  fourth  order  Runge-Kutta  integration  is  normally  used,  a  variable 
order  and  step  size  predictor-corrector  method  (Ref.  14)  was  also  used,  when  checking  integra¬ 
tion  accuracy. 


3.1  Effect  of  Number  of  ErscMises  Matched — No  Noise 

In  order  to  assess  the  >-elative  importance  of  each  of  the  measured  time  histories,  various 
combinations  of  two  and  three  elements  as  well  as  the  full  four  elements  of  the  output  vector 
(equation  (15))  were  matched.  In  the  absence  of  measurement  noise  the  weightings  given  to  the 
records  (i.e.  the  diagonal  elements  of  R)  were  fixed  approximately  inversely  proportional  to  the 
square  of  a  typical  expected  noise  level.  The  results  are  in  fact  fairly  insensitive  to  variations  in 
relative  weightings. 

In  order  to  avoid  the  correlation  problem  discussed  in  Section  2.4  the  biases  bv,  bx  and  b„ 
were  fixed  at  their  correct  values  thus  giving  a  six-element  parameter  vector  (equation  (16)). 
The  manoeuvre  shape  used  (Manoeuvre  1)  is  produced  by  the  inputs  shown  in  Figure  2.  Starting 
from  steady  level  flight  the  manoeuvre  assumes  a  linear  variation  of  longitudinal  acceleration, 
ax,  and  periodic  variations  of  normal  acceleration,  az,  and  pitch  rate,  q.  The  output  airspeed, 
incidence  and  pitch  attitude  over  a  60-second  time  span  are  shown  in  Figure  3,  although  only 
the  first  20  seconds,  at  a  sampling  rate  of  10  per  second,  was  used  at  this  stage.  Figure  3  shows 
both  the  "measured”  records  and  those  calculated  assuming  zero  biases  but  with  initial  values 
matched. 

A  summary  of  results,  each  obtained  after  four  iterations  and  a  typical  processing  time  of 
60  seconds  on  the  PDPIO  computer,  is  presented  in  Table  I.  An  indication  of  the  relative 
accuracy  of  the  estimates  is  provided  by  the  number  in  parentheses  which  are  the  “Cramer-Rao" 
bounds  computed  with  the  fixed  R  matrix.  It  is  clear  from  the  results  that  the  use  of  three  mat  ;hed 
time  histories  can  produce  almost  as  good  results  as  four  matched  records.  This  is  particularly 
true  with  the  V ,  <*,  6  records  which  produce  results  for  the  estimated  parameters,  their  confidence 
bounds  and  time  history  match  errors  (expressed  as  RMS  error)  in  close  agreement  with  the  full 
match  case.  The  use  of  only  two  time  histories  leads  to  a  general  deterioration  of  results  with  the 
V,  ot  and  V,  6  cases  being  the  least  unacceptable. 

It  is  concluded  from  this  study  that  matching  the  height,  h,  record  only  provides  a  marginal 
improvement  of  the  results  compared  to  a  V,  a,  6  match.  Considering  also  the  difficulties  of 
obtaining  accurate  height  measurements  and  the  favourable  reduction  in  the  system  equations 
when  li  is  neglected  (Section  2.3),  height  has  been  removed  from  the  system  equations  in  the 
following  section*. 


3.2  Effect  of  Parameter  Correlations — No  Noise 

In  this  section  all  nine  unknown  parameters  are  estimated  (equation  (16))  with  measure¬ 
ments  assumed  to  be  noise  free.  Because  of  the  relations  between  the  initial  conditions  and  several 


Number  of  Time  Histories  Matched — No  Noise 


bv,  ba,  be  fixed  at  correct  value. 

Four  iterations,  typical  CPU  time  about  60  s. 

A  priori  parameter  vector  =  [0  0, 0  0, 0  0,  101  0,  10-3,  0  01]T. 
Manoeuvre  1. 


of  the  biases,  a  number  of  different  procedures  was  tried  as  discussed  in  Section  2.4  and 
Appendix  3.  Six  methods  used  are  summarised  in  the  table  below. 


Method  1  uses  the  unmodified  sensitivities  summarised  in  Appendix  1  while  Methods  2-5 
include  modifications  described  in  the  corresponding  sections  of  Appendix  3.  Method  6  proceeds 
as  in  Method  1  but  leaving  out  the  initial  conditions  from  the  parameter  vector.  However, 
initial  conditions  are  adjusted  at  each  iteration  using  the  relations  in  Section  2.4.  Thus  Method  6 
is  like  Method  4  but  without  the  sensitivity  matrix  elements  being  modified.  This  appears  to  be 
the  approach  used  in  Reference  9. 

The  manoeuvre  used  is  the  same  as  in  the  previous  section  (Manoeuvre  1)  as  are  the  fixed 
elements  of  the  weighting  matrix.  Twenty  seconds  of  record  at  20  samples  per  second  was 
analysed  and  typical  computing  time  was  400  seconds  for  20  iterations.  For  this  and  all  subse¬ 
quent  computations  the  a  priori  values  for  the  biases  were  set  to  zero  and  those  for  the  initial 
conditions  u( 0),  n(0)  and  0(0)  set  to  values  determined  from  the  “measurements”. 

A  comparison  of  the  performance  of  the  six  methods  with  noise  free  data  is  shown  in 
Table  2.  With  Method  1,  although  the  “Cramer-Rao”  bounds  are  relatively  small  the  estimated 
parameter  values  are  clearly  biased.  Examination  of  the  normalised  correlation  matrix  reveals 
a  number  of  highly  correlated  parameters,  as  expected,  the  worst  being  correlations  between 
bx-w( O)-0(O)  and  between  bv-u{0).  An  immediate  improvement  in  parameter  values  is  obtained 
with  Method  2  when  0(0)  is  removed  from  the  list  of  independent  parameters.  Further  improve¬ 
ments  are  achieved  with  Methods  3  and  4  which  give  very  similar  results,  as  may  be  expected, 
but  Method  5  shows  little  consistent  improvement  over  the  previous  two  methods.  While 
estimated  parameter  values  improve  going  from  Method  1  to  Method  5  it  is  worth  noting  that 
the  confidence  bounds  for  Methods  2-5  are  worse  than  that  of  Method  1.  Finally  Method  6  is 
a  poor  performer  giving  strongly  biased  results.  The  “Cramer-Rao”  bounds  appear  to  be  very 
small  in  this  case  but  it  should  be  remembered  that  the  problem  formulation  is  inconsistent. 


3.3  Effect  of  Noise,  Record  Length  and  Sampling  Rate 

The  effects  of  both  measurement  and  process  noise  on  the  estimation  procedures  of  the 
previous  section  were  investigated  using  different  (Gaussian)  noise  levels.  Table  3  below  lists 
the  noise  standard  deviations  on  each  of  the  measurement  channels  for  five  different  noise 
levels  used. 

The  first  two  columns  in  Table  3  are  measurement  noise  only  cases,  the  first  column  repre¬ 
senting  a  fairly  high  noise  level  and  the  second  a  low  noise  level  similar  to  that  used  in  Reference  9. 
The  last  three  columns  include  both  process  and  measurement  noise.  Level  3  is  a  low  process 
noise  case  (Ref.  9)  confined  with  low  measurement  noise;  Level  4  combines  low  process  noise 
with  high  measurement  noise  and  Level  5  has  high  noise  levels  in  both  process  and  measurement 
noise.  Even  the  low  noise  levels  in  Table  3  are  generally  considerably  higher  than  the  levels  des¬ 
cribed  in  Reference  10,  which  correspond  to  highly  accurate  instrumentation.  Instrumentation 
of  such  accuracy  is  not  yet  in  common  use  for  flight  testing. 
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N  =  400,  record  length  =  20  s,  R  =  diag[l  0,  250000,  750000]fixed. 
Only  K,  a,  0  time  histories  matched 
20  iterations,  typical  CPU  time  about  400  s. 

A  priori  parameter  vector  =  [0  0, 0  0, 0  0, 0  0, 0  0, 0  0,  99  431,  17 
Manoeuvre  1. 


TABLE  3 
Noise  Levels 


Standard 

error 

Measurement  noise  only 

Process  and  measurement  noise 

Level  1 

Level  2 

Level  3 

Level  4 

Level  5 

a(a.v),  m/s2 

0  0 

00 

f 

<r(ar),  m/s2 

0  0 

0  0 

BUS 

0  1 

a(q),  rad/s 

00 

0  0 

O' 001 

0  001 

<j(T),  m/s 

10 

0  1 

01 

10 

cr (*),  rad 

0  002 

0  001 

0  001 

0  002 

o(0),  rad 

0  01 

0  001 

0  001 

0  01 

The  effects  of  measurement  noise  on  the  results  of  the  previous  section  are  summarised  in 
Table  4.  The  test  details  remain  the  same  except  that  Level  I  noise  has  been  added  to  the  measure¬ 
ments  and  consequently  the  estimate  for  the  measurement  noise  matrix,  R,  is  updated  in  the 
course  of  the  analysis  as  explained  in  Section  2.1.  Thus  the  Cramer-Rao  bound  (equation  10) 
should  be  a  good  estimate  of  the  accuracy  of  the  estimates.  The  results  in  Table  4  represent  the 
average  values  plus  or  minus  the  standard  deviation  (or  scatter)  obtained  from  8  to  10  separate 
runs.  The  scatter  should  be  comparable  with  the  Cramer-Rao  bounds  shown  in  parentheses. 
In  comparing  Tables  2  and  4,  the  results  for  Method  1  (all  parameters  identified)  are  seen  to  be 
very  similar.  The  scatter  is  about  twice  the  Cramer-Rao  bound  for  most  parameters,  indicating 
room  for  improvement.  However,  the  mean  values  for  all  the  parameters  are  within  about 
one  Cramer-Rao  bound  of  their  true  value.  The  results  from  Method  6  remain  erratic  as  before 
although  apparently  providing  a  reasonable  fit.  The  Cramer-Rao  bounds  in  this  case  give  little 
indication  of  the  accuracy  of  the  estimated  parameter  values.  The  remaining  methods  all  show 
a  marked  deterioration  in  the  presence  of  measurement  noise.  The  average  results  no  longer 
compare  favourably  with  those  from  Method  I  and  the  scatter  bands  are  considerably  worse. 
Method  4,  not  shown  in  Table  4,  provides  results  very  similar  to  Method  3.  The  deterioration  in 
performance  of  Methods  2  to  5  in  the  presence  of  noise  is  probably  linked  to  the  fact,  as  discussed 
in  Section  2.4,  that  the  relations  assumed  by  equation  (23)  no  longer  hold  exactly. 

It  is  concluded  that,  despite  the  correlations  between  several  of  the  parameters.  Method  1 
is  likely  to  provide  the  best  results.  In  order  to  optimise  further  the  performance  of  Method  I 
the  effects  of  record  length,  sampling  time  and  measurement  noise  level  have  been  examined 
and  the  results  are  shown  in  Table  5.  These  are  average  results  over  7  to  10  simulations  with 
the  scatter  and  Cramer-Rao  bound  (in  parentheses)  also  shown.  Record  lengths  vary  from 
20  to  60  seconds  and  the  number  of  samples,  N,  corresponds  to  sampling  rates  of  either  20  or 
40  samples  per  second.  The  first  five  columns  of  results  are  for  Level  1  (high)  measurement  noise, 
the  first  (reference)  column  being  a  repeat  of  the  Method  l  results  from  Table  4.  Column  2  is 
for  40  seconds  of  record  at  the  same  sampling  rate  and  shows  a  notable  improvement  in  estimated 
parameters.  The  scatter  of  the  results  is  much  reduced  and  shows  good  agreement  with  the 
Cramer-Rao  bounds  in  all  cases.  However,  the  scatter  in  some  cases,  in  particular  />„  is  still 
somewhat  large.  Doubling  the  sampling  rate  (column  3)  produces  further  improvement  in  both 
estimates  and  scatter  as  does  an  increase  in  the  record  length  from  40  to  60  seconds  (column  4) 
at  the  lower  sar  p'ine  rate.  Column  5  gives  results  for  60  seconds  of  record  at  the  higher  (40  per 
second)  sampling  rates,  i  he  improving  results  with  increasing  record  length  are  partly  due  to 
decreasing  correlations  between  parameters,  but  close  examination  of  the  results  suggests  that 
there  may  still  be  some  bias  in  bax  and  perhaps  also  0(0). 

The  last  two  columns  of  Table  5  are  for  lower  (level  2)  measurement  noise.  As  expected 
results  are  improved  for  both  estimated  values  and  scatter.  At  this  noise  level  40  seconds  of 
record  is  sufficient  to  identify  all  the  biases  to  within  about  IOn0  or  better. 

The  addition  of  process  noise  can  be  expected  to  affect  the  results  adversely,  since  it  has  not 
been  accounted  for,  and  some  results  are  shown  in  Table  6.  These  can  be  compared  directly  with 
the  relevant  columns  of  Table  5.  For  example  columns  I  and  2  of  Table  6  should  be  compared 
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Results  are  averages  of  7-10  runs. 
Manoeuvre  1. 


with  columns  7  and  5,  respectively,  of  Table  5.  While  the  respective  Cramer-Rao  bounds  remain 
the  same,  the  scatter  increases  markedly.  This  deterioration  is  relatively  much  worse  in  the  case 
with  low  measurement  noise  (column  1).  Nevertheless  most  of  the  biases  are  reasonably  well 
identified  when  60  seconds  of  record  are  used,  the  main  exception  being  ba  which  has  a  large 
scatter  band.  When  either  the  record  length  is  decreased  (column  3)  or  the  level  of  process  noise 
is  increased  (column  4)  the  accuracy  of  estimation  of  several  of  the  other  parameters  also  becomes 
unacceptable,  but  it  is  clear  that  bq  is  always  well  identified  followed  probably  by  baz- 

TABLE  6 

Effect  of  Process  and  Measurement  Noise — Method  1 


Para¬ 

meter 

True 

value 

bax,  m/s2 

0-1 

bai ,  m/s2 

0-1 

bq,  rad/s 

0-002 

bv,  m/s 

1-0 

ba,  rad 

0-002 

be,  rad 

0  01 

w(0),  m/s 

98-48 

K-(O),  m/s 

17-36 

0(0),  rad 

0  175 

Time,  sec 
N 

Noise 

■ 

0-086±0  006 

(0  0020) 

0 -099  ±0  005 
(0  0005) 

0-0020  ±0-00002 
(0-000001) 
1-038  ±0134 
(0-025) 

0 -0019±0-0015 

(0  0002) 

0  0114±0  0004 
(0-0002) 
98-436±0- 172 
(0-026) 
17-360±0- 161 
(0-020) 

0- 173±0  001 
(0-0002) 

60 
2400 
Level  3 


2 

3 

0  087  ±0-011 

0-051  ±0-020 

(0-0051) 

(0-0125) 

0-095±0-003 

0-092  ±0-01 1 

(0-0019) 

(0  0039) 

0-0020 ±0-00002 

s 

6 

-H 

8 

8 

6 

(0-00001) 

(0  00002) 

1  •  132±0- 163 

1  213±0  526 

(0-135) 

(0-230) 

0-0018±0-0026 

0  0052  ±0  0029 

(0-0006) 

(0-0013) 

0-01 15±0-001 1 

0-0158  ±0-0023 

(0-0006) 

(0-0013) 

98  *  3 1 5  ±0  - 197 

98 -303  ±0-520 

(0-147) 

(0-238) 

17-379±0-306 

17-028±0-340 

(0-073) 

(0-140) 

0-1728  ±0  001 

0-1690  ±0  002 

(0-0005) 

(0  0013) 

60 

40 

2400 

1600 

Level  4 

Level  4 

0  089  ±0-02l 
(0  0053) 
0-092  ±0-008 
(0-0020) 

0-0020 ±0-00005 
(0-00001) 
1-234  ±0-479 
(0-138) 

0  •  000 1 7  ±0  0034 
(0-0006) 

0-01 11  ±0-0021 
(0  0006) 
98-219±0-556 
(0-150) 

17 -505  ±0-354 
(0  075) 

0  1731  ±0  002 
(00005) 

60 
2400 
Level  5 


Only  V,  a ,  8  records  matched;  R  =  diag[l  0,  250000,  10000]  initially. 

10  iterations — CPU  time  about  20  min  for  cases  1,  2,  4  and  14  min  for  case  3. 

A  priori  parameter  vector  =  (0-0,  0-0, 0-0,  0-0,  0-0,  0-0,  99-431,  17-734,  0-  185]T. 
Results  averages  of  seven  runs. 

Manoeuvre  1. 


3.4  Effect  of  Manoeuvre  Shape 

All  of  the  results  so  far  presented  were  obtained  using  the  manoeuvre  shape  (Manoeuvre  1) 
shown  in  Figure  2.  It  is  known  that  manoeuvre  shape  influences  the  amount  of  information  in 
the  measurements  and  efforts  have  been  made  elsewhere  in  the  development  of  optimal  manoeuvre 
shapes  (see  e.g.  Ref.  15).  No  attempt  is  made  here  to  develop  an  optimal  shape  but  in  this  Section 
results  with  changing  manoeuvre  shapes  are  reported  so  as  to  bring  out  the  effect  this  can  have  on 
the  results. 

The  shapes  used  are  as  follows: 

Manoeuvre  1. — Shown  in  Figure  2,  the  inputs  consist  of  a  longitudinal  deceleration  at 
0- 1  m/s2,  an  oscillatory  normal  acceleration  with  an  amplitude  of  2  m/s2  and  an  oscillatory  pitch 
rate  with  an  amplitude  of  0-2  rad.  The  frequency  of  oscillation  in  each  case  is  2  rad/s. 


Effect  of  Manoeuvre  Shape-Method  1 
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/f  prior/  values,  biases  =0  0,  initial  conditions  from  measurements. 
Results  with  process  noise  are  average  of  nine  runs. 


Manoeuvre  2. — Shown  in  Figure  4,  the  only  difference  from  the  previous  inputs  being  in 
the  normal  acceleration  which  has  now  been  made  symmetrical  about  the  undisturbed  value. 
The  output  airspeed,  incidence  and  attitude  are  shown  in  Figure  5.  Level  1  (measurement)  noise 
is  shown  on  the  outputs  and  the  calculations  assume  zero  biases.  Thus  the  differences  between 
the  two  sets  of  curves  represent  the  effects  of  the  biases.  The  major  difference  between  this  and 
the  previous  manoeuvre  is  in  the  angle  of  attack  which  no  longer  reaches  the  very  high  values 
caused  by  the  asymmetric  normal  acceleration  input. 

Manoeuvre  3. — This  manoeuvre  is  very  similar  to  the  previous  except  that  the  normal 
acceleration  amplitude  has  been  increased  to  10  m/s2  thus  providing  a  wider  range  in  angle 
of  attack. 

Manoeuvre  4  (Fig.  6). — This  manoeuvre  differs  from  the  previous  only  in  having  a  different 
oscillation  frequency  for  the  input  normal  acceleration.  The  differences  in  the  output  records 
are  shown  in  Figure  7. 

Manoeuvre  5  (Fig.  8). — This  is  a  simulation  of  a  real  manoeuvre  produced  by  a  pilot  applying 
stick  inputs  to  produce  a  roller  coaster  type  of  motion,  i.e.  push-over/pull-up/push-over.  The 
resulting  accelerations  and  pitch  rate  are  similar  to  those  shown  in  Figure  8,  with  the  higher 
frequency  oscillations,  particularly  in  pitch  but  also  in  normal  acceleration,  being  due  to  the 
Short  Period  response.  The  outputs,  both  measured  (with  Level  1  noise)  and  calculated  (assum¬ 
ing  zero  biases)  are  shown  in  Figure  9. 

The  results  using  these  five  manoeuvre  shapes  are  shown  in  Table  7.  For  this  study  40 
seconds  of  record  at  40  samples  per  second  were  used.  Consider  first  the  five  columns  of  results 
with  Level  1  measurement  noise  only.  For  Manoeuvre  1  the  results  are  a  repeat  of  those  in 
Table  5,  column  3.  With  Manoeuvre  2  the  smaller  angle  of  attack  range  leads  to  a  general 
worsening  of  results,  both  in  parameter  value  and  Cramer-Rao  bound.  Of  the  biases,  only  bq 
and  perhaps  baz  can  be  taken  as  satisfactorily  identified.  With  Manoeuvre  3  the  larger  disturbance 
amplitude  produces  a  general  improvement  in  results  all  round.  However,  bg  and  bax  still  are 
strongly  biased  as  indicated  by  the  Cramer-Rao  bounds.  Manoeuvre  4  introduces  two  dominant 
frequencies,  particularly  obvious  in  the  angle  of  attack  output  (Fig.  7\and  this  extra  information 
leads  to  a  further  general  improvement  in  the  parameter  estimates  and  error  bounds.  In  addition 
biases  in  bax  and  be  are  no  longer  obvious.  Finally,  the  information  content  of  Manoeuvre  5  is 
such  as  to  produce  the  best  results  in  this  Section.  Only  bx  (and  the  related  initial  value,  iv(0)) 
still  appear  troublesome. 

The  last  two  columns  of  Table  7  include  both  process  and  measurement  noise  with 
Manoeuvre  5.  The  Level  5  (high  process  and  measurement  noise)  results  are,  on  average,  very 
similar  to  the  measurement  noise  only  (Level  I)  case  although  the  scatter,  based  on  the  standard 
deviation  over  nine  runs,  is  two  to  three  times  worse  than  the  Cramer-Rao  bound  now.  The 
Level  3  (low  process  and  measurement  noise)  results  produce  a  marked  improvement  with  all 
parameters  being  well  identified.  Although  the  scatter  is  about  five  times  the  Cramer-Rao 
bound,  it  is  comparable  and  in  some  cases  better  than  the  Cramer-Rao  bounds  obtained  with 
the  higher  levels  of  measurement  noise.  In  particular  />a  and  >r(0)  are  considerably  improved. 

3.5  Summary  of  Results 

The  results  for  the  longitudinal  kinematic  system  studied  here  have  suggested  that  the  use 
of  the  three  elements  V ',  a,  0  of  the  output  vector  produce  results  as  good  as  those  using  the  full 
vector,  including  height.  Reduced  computation  time  and  simplicity  also  favour  a  reduced  system. 

The  problem  of  correlations  between  the  state  initial  conditions  and  some  parameters  was 
investigated  using  several  possible  methods  proposed  to  account  for  the  correlations.  Although 
improved  performance  was  obtained  from  these  methods  in  the  noise  free  case,  the  addition  of 
measurement  noise  led  to  marked  deterioration.  Best  results,  with  measurement  noise,  were 
obtained  when  all  parameters  including  initial  conditions  were  independently  extracted. 

Correlations  decreased  as  record  lengths  increased  with  consequent  decrease  in  the  bias  in 
some  estimates  and  decrease  in  the  scatter  band.  Reasonable  accuracy  could  be  obtained  with  a 
record  length  of  40  seconds  at  40  samples  per  second.  Reduction  of  measurement  noise  levels 
was  also  found  to  lead  to  marked  improvement  in  parameter  estimates.  The  main  effect  of  process 
noise  was  an  increase  in  the  scatter  of  results.  Also  bias  in  some  estimates  appeared  where  previ¬ 
ously  absent  with  measurement  noise  alone. 
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It  was  clear  that  the  manoeuvre  shape  had  an  important  influence  on  the  results.  A  practical 
push-over/pull-up  manoeuvre  was  found  to  provide  satisfactory  results  for  most  parameters 
with  a  record  length  of  40  seconds,  particularly  for  lower  noise  levels  (both  measurement  and 
process).  Further  improvements  may  well  be  attainable  by  developing  a  manoeuvre  specially 
suited  to  the  purpose  of  identifying  instrument  biases. 


4.  CONCLUSIONS 

A  general  Maximum  Likelihood  program,  suitable  for  estimation  of  parameters  in  non¬ 
linear  systems,  has  been  described.  For  successful  application  to  a  particular  system,  process 
noise  levels  should  be  small.  A  summary  of  the  theoretical  background  has  been  given  and  the 
program  structure  has  been  developed  via  subroutines  to  allow  ready  modification  to  suit  the 
problem  in  hand. 

The  program  has  been  applied  here  to  a  study  involving  longitudinal  aircraft  kinematics. 
In  this  case  the  kinematic  equations,  treated  separately  from  a  dynamic  representation  of  air¬ 
craft  motions,  are  free  from  process  noise  to  a  good  approximation.  This  separate  treatment 
of  the  kinematics  can  be  used  for  checking  the  compatibility  of  measurements  and  for  estimating 
instrument  systematic  errors.  Such  an  approach  is  useful  for  reconstructing  “error-free”  records 
prior  to  dynamic  analysis  aimed  at  estimating  aerodynamic  parameters,  and  is  particularly 
valuable  when  the  dynamic  model  is  uncertain  and/or  non-linear. 

The  study,  using  simulated  data,  has  aimed  at  establishing  the  conditions  under  which  the 
Maximum  Likelihood  algorithm  can  successfully  extract  the  instrument  biases,  and  to  set  a 
baseline  for  comparison  with  other  methods.  It  has  been  found  that  the  present  method  is  feasible 
provided  record  lengths  and  sampling  rates  are  adequate.  The  importance  of  using  high  quality 
instrumentation  to  minimise  noise  has  also  been  brought  out.  Further,  it  has  been  shown  that 
the  manoeuvre  shape  can  have  a  considerable  influence  on  the  quality  of  results  and  an  effort 
to  develop  an  optimal  manoeuvre  shape  may  be  worthwhile.  For  the  present  system,  using  a 
practical  manoeuvre  shape,  40  seconds  of  record  at  40  samples  a  second  with  reasonably  low 
measurement  and  process  noise  levels  was  successful  in  extracting  most  parameters  to  within 
10%  or  better  of  their  true  values.  The  precise  levels  of  accuracy  acceptable  would  depend,  to 
some  extent,  on  the  requirements  of  the  subsequent  analysis  and  would  need  to  be  investigated 
further.  Finally,  the  performance  needs  to  be  checked  with  real,  as  opposed  to  simulated,  flight 
test  data. 
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APPENDIX  1 

Sensitivity  Matrix  (Subroutine  SENS) 


The  parameter  vector  has  dimension  9: 

f  =  [baz,  baz ,  bt,  bv,  ba,  be,  u{ 0),  h(0),  6(0)]t 
and  the  output  vector  has  dimension  4: 

=  [I^outj  aVouti  60ut,h0ut]T. 

Using  the  output  definition  in  equation  (15),  the  transpose  of  the  4x9  sensitivity  matrix 
becomes : 
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where 

Vout.lj  =  (“  •  u( j  +  W  .  W(j)l(u*+  H’2)1'2, 

«vou< ,ij  =  («  •  wij-[w-(<frn+bQ)xa]uu)l(u2+[w-(qm+bg)xa]2), 
6out.(j  =  0[j, 
flout. £,  —  fl(j- 


The  above  relations  hold  for  those  elements  of  the  parameter  vector  as  listed  above  in 
the  sensitivity  matrix.  The  only  exception  is 

“voul,6,  =  (u  .  wbq-[w-(qm+be)xa]ub<l-u .  xa)l(u2+[^~(<]m+bq)xa]2). 


APPENDIX  2 


Differential  Equations  (Subroutine  DERIVS) 

In  order  to  calculate  the  state  (equation  (14))  and  the  elements  of  the  sensitivity  matrix 
(Appendix  I)  the  following  system  of  ordinary  differential  equation  needs  to  be  integrated  with 
respect  to  time: 

u  =  —  (qm+bq)n +(axm-\-bax)~- g  sin  # 
w  =  (qm+bq)u+(a:m+ba2)+gcos  0 
&  =  <7m  ~rbq 
h  =  u  sin#  —  tf  cos  6 

Ubajc  =  (qn\\bq)u'tiay  +  I 
'''fax  ~  (#m  f  bq)Ut,llx 

hbax  =  sin# .  Ubax— COS  #  .  Ubujt 

Ufa;  ~  ~(qm-\-bq)Wba: 

''fa?  =  (9mJrbq)Uba;Jr  1 

llbj;  =  sin  #  .  Uba,  COS  #  .  H(,a. 

%  =  ~(qmJrbq)\\bq  —  ir — g  cos  #.  0bq 
Wbq  =  (qm +bq)ubq-ru-g  sin  #  .  9bq 
6hq  =  1  0 

hbq  =  sin  #  .  ubq— cos  6  .  \\  bq+(u  cos  #+tr  sin  6)0bq 
#«<<*>  =  (qm-]~bq)\\'u{0) 

'Vu(O)  =  (^tn  +  6<j)Uii(0) 
b mo  =  sin  # .  mH(0)  — cos  #  .  tt'u (0) 

#«(0|  =  — (<7m  + 
wm(0)  =  (qm-\-bq)uu <0) 

Ufr(O)  =  sin  #  Uu'(O)  cos  #  . 

#s<0|  —  (<7m  +  (fy)>‘'0<(i)  — g  COS  # 

w>(0,  =  (^m  +  ^M^Oi-g  sin  # 

/>»(0i  =  sin  #  .  M#(0|  — cos  #  .  11-0,01  +  h  cos  #+»»•  sin  # 


APPENDIX  3 


Effect  of  Correlations  on  Sensitivity  Matrix 

As  explained  in  Section  2.4,  when  relationships  between  initial  conditions  and  biases  are 
accounted  for,  it  is  possible  to  reduce  the  number  of  parameters  to  be  estimated.  At  the  same 
time,  alterations  need  to  be  made  to  some  elements  of  the  sensitivity  matrix  defined  in  Appendix  I . 
Four  possible  approaches  which  have  been  examined  and  the  resulting  changed  elements  of  the 
sensitivity  matrix  are  summarised  below. 

(1)  Parameter  vector,  [bax,baz,bq,bv,bx,u(0),w(0),d(0)]T — bg  treated  as  a  function  of  0(0). 
Changed  sensitivity  elements: 

0out.fl(O)  =  0 

(2)  Parameter  vector,  l6ax,£aji6<,,w(O),u(O),0(O)]T  —  bv,ba,bg  treated  as  functions  of  w(0), 
h(0),  0(0). 

Changed  sensitivity  elements: 

F out.utO)  =  (u  •  Wu(0|  +  M’  •  H|1(o,)(/m2  +  M'2)1  2— m(0)/(m2(0)  +  H-2(0))I/2 

Fout.«(0l  =  («  ■  Wu>(0|  +  W  .  «',c(0|)/(M2  +  U’2)1  2  —  »*’(0)/«2(0)+  M'2(0))'' 2 
avOU|,U(0)  =  (u  .  H«,0)-[u'-(^m+^).V1]ttM,0))/(U2+[H’-(^m  +  N)jtc,]2) 

+(H(0)  -  (qm  +  bq).\x)l(u-(0)  r  [h(0)  -  (qm( 0) + bq)xxf) 

avou,.i<(0)  =  (w.  >««-(0)-[>‘  -(?m  +  /»<f).V;I]MK(0.)/(M2+[M—  (<?m+6?)jt«]2) 

—u(0)/(w2(0)  +  (h  (0)  ~  (qm(0)+bq)x^) 

0ou(.tf(O)  =  0 

(3)  Parameter  vector,  [baz,baz,bq,bv,bx,bg]T  —  w(O),n(O),0(O)  treated  as  functions  of  bv ,  bx,  bg 
Changed  sensitivity  elements: 

Fout.fty  =  1  --[(«  ■  Mu(0|  +  U-  .  H'U(0))COS  a(0) 

+  (u  •  w«(0)4-t>' .  u  «.(o,)sin  a(0)]/(w2  +  w2)1  2 
F out,&a  =  [(w  ■  Wk(0)  "T t*‘ .  >t  „(oi)sin  *(0) 

-(«  rtf  .  M„.(o,)cos  a(0)] .  (u2(0)  +  »r2(0))1  2/(u2-f-»r2)1  2 

F i out.fts  =  — (u  .  •  "  @(oi)/(u2-t- ii-)1  - 

avoul.hv  =  [(»'  •  .  H  a(0,)cos  a(0)  f(u- .  .  wy,0»)sin  a(0)]/(u2+ir2) 

avoul,&i  =  1  [(w  .  u‘u(0i  ii‘  ■  w«<oi)sin  3^(0) 

-r (»'  •  w«<0)  -u  .  H-,,|0i)cos  *(0)] .  (w2(0) -t- •i'2(0))1  2/(m2  +  h'2) 
avoul.*fl  =  •  M«<0(  M  .  »r#,o,)/(H-  —I*'-) 

0out ,bo  —  0 

These  relations  assume  initially  steady  flight,  i.e.  pitch  rate  equals  zero. 


(4)  Parameter  vector,  [bat,bt,u(0),w(0),d(0)]—bv,ba,be  and  bax  treated  as  functions  of  u(0), 
h(0),  fl(0). 

Changed  sensitivity  elements : 

^out.uiOi,  l7 out,u'(0)>  avou,,u(0),  avout,uj(0)  and  ^out.siO)  same  as  listed  in  approach  2.  In 
addition: 

^01 K,tf(0)  =  [u.  uet0)+w .  wet0)  — (azm(0)+ioz)sec2  0(O)(m.  Ubax+w.  Wbax)V(u2jrw2y  2 
“voul.eiO)  =  [u  ■  -  Mtf(0)-(azm(0)  +  iaj)sec2  0(O)(u  .  H'bax  —  H’  .  tO]/(“2  +  »*2) 

Initially  steady  flight,  i.e.  zero  pitch  rate,  has  been  assumed  in  these  relations. 
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aircraft  dynamic  flight  test  data  has  been  studied.  Using  simulated  time  histories  of  longitudinal 
manoeuvers,  the  conditions  for  satisfactory  performance  have  been  identified.  It  has  been  shown 
that  for  a  practical  maneuver  shape,  record  length  and  sampling  rate  and  for  reasonable  noise 
levels  on  the  measured  data,  instrument  errors  can  be  identified  to  good  accuracy  and  a  set  of 
compatible  time  histories  reconstructed. 
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